home *** CD-ROM | disk | FTP | other *** search
/ Inside Multimedia 1994 April / Inside Multimedia CD-ROM (April 1994).iso / prg / gs / gssource.exe / GXPCOPY.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-08-24  |  16.9 KB  |  504 lines
  1. /* Copyright (C) 1992 Aladdin Enterprises.  All rights reserved.
  2.    Distributed by Free Software Foundation, Inc.
  3.  
  4. This file is part of Ghostscript.
  5.  
  6. Ghostscript is distributed in the hope that it will be useful, but
  7. WITHOUT ANY WARRANTY.  No author or distributor accepts responsibility
  8. to anyone for the consequences of using it or for whether it serves any
  9. particular purpose or works at all, unless he says so in writing.  Refer
  10. to the Ghostscript General Public License for full details.
  11.  
  12. Everyone is granted permission to copy, modify and redistribute
  13. Ghostscript, but only under the conditions described in the Ghostscript
  14. General Public License.  A copy of this license is supposed to have been
  15. given to you along with Ghostscript so you can know your rights and
  16. responsibilities.  It should be in a file named COPYING.  Among other
  17. things, the copyright notice and this notice must be preserved on all
  18. copies.  */
  19.  
  20. /* gxpcopy.c */
  21. /* Path copying and flattening for Ghostscript library */
  22. #include "math_.h"
  23. #include "gx.h"
  24. #include "gserrors.h"
  25. #include "gxfixed.h"
  26. #include "gxarith.h"
  27. #include "gzpath.h"
  28.  
  29. /* Forward declarations */
  30. private int copy_path(P4(const gx_path *, gx_path *, fixed, int));
  31. private int flatten_recur(P8(gx_path *,
  32.   fixed, fixed, fixed, fixed, fixed, fixed, fixed));
  33. private int flatten_sample(P8(gx_path *, int,
  34.   fixed, fixed, fixed, fixed, fixed, fixed));
  35.  
  36. /* Copy a path */
  37. int
  38. gx_path_copy(const gx_path *ppath_old, gx_path *ppath, int init)
  39. {    return copy_path(ppath_old, ppath, max_fixed, init);
  40. }
  41.  
  42. /* Flatten a path. */
  43. /* If flatness is zero, use sampling rather than subdivision: */
  44. /* this is important for Type 1 fonts. */
  45. private fixed
  46. scale_flatness(floatp flatness)
  47. {    fixed scaled_flat = float2fixed(flatness);
  48.     return (scaled_flat > int2fixed(100) ? int2fixed(100) :
  49.         scaled_flat < fixed_half ? fixed_0 :
  50.         scaled_flat);
  51. }
  52. int
  53. gx_path_flatten(const gx_path *ppath_old, gx_path *ppath, floatp flatness, int in_BuildChar)
  54. {    return copy_path(ppath_old, ppath,
  55.              (in_BuildChar ? fixed_0 : scale_flatness(flatness)),
  56.              1);
  57. }
  58.  
  59. /* Add a flattened curve to a path. */
  60. int
  61. gx_path_add_flattened_curve(gx_path *ppath,
  62.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3,
  63.   floatp flatness)
  64. {    return flatten_recur(ppath, x1, y1, x2, y2, x3, y3,
  65.                  scale_flatness(flatness));
  66. }
  67.  
  68. /* Copy a path, optionally flattening it. */
  69. /* If the copy fails, free the new path. */
  70. private int
  71. copy_path(const gx_path *ppath_old, gx_path *ppath, fixed scaled_flat,
  72.   int init)
  73. {    gx_path old;
  74.     const segment *pseg;
  75.     int code;
  76. #ifdef DEBUG
  77. if ( gs_debug['p'] )
  78.     gx_dump_path(ppath_old, "before copy_path");
  79. #endif
  80.     old = *ppath_old;
  81.     if ( init )
  82.         gx_path_init(ppath, &ppath_old->memory_procs);
  83.     pseg = (const segment *)(old.first_subpath);
  84.     while ( pseg )
  85.        {    switch ( pseg->type )
  86.            {
  87.         case s_start:
  88.             code = gx_path_add_point(ppath, pseg->pt.x, pseg->pt.y);
  89.             break;
  90.         case s_curve:
  91.            {    curve_segment *pc = (curve_segment *)pseg;
  92.             if ( scaled_flat == max_fixed )    /* don't flatten */
  93.                 code = gx_path_add_curve(ppath,
  94.                     pc->p1.x, pc->p1.y,
  95.                     pc->p2.x, pc->p2.y,
  96.                     pc->pt.x, pc->pt.y);
  97.             else
  98.                 code = flatten_recur(ppath,
  99.                     pc->p1.x, pc->p1.y,
  100.                     pc->p2.x, pc->p2.y,
  101.                     pc->pt.x, pc->pt.y,
  102.                     scaled_flat);
  103.             break;
  104.            }
  105.         case s_line:
  106.             code = gx_path_add_line(ppath, pseg->pt.x, pseg->pt.y);
  107.             break;
  108.         case s_line_close:
  109.             code = gx_path_close_subpath(ppath);
  110.             break;
  111.            }
  112.         if ( code )
  113.            {    gx_path_release(ppath);
  114.             if ( ppath == ppath_old )
  115.                 *ppath = old;
  116.             return code;
  117.            }
  118.         pseg = pseg->next;
  119.     }
  120.     ppath->position = old.position;        /* restore current point */
  121. #ifdef DEBUG
  122. if ( gs_debug['p'] )
  123.     gx_dump_path(ppath, "after copy_path");
  124. #endif
  125.     return 0;
  126. }
  127. /* Internal routine to flatten a curve. */
  128. /* This calls itself recursively, using binary subdivision, */
  129. /* until the approximation is good enough to satisfy the */
  130. /* flatness requirement.  The starting point is ppath->position, */
  131. /* which gets updated as line segments are added. */
  132.  
  133. /* Maximum number of points for sampling if we want accurate rasterizing. */
  134. /* (num_sample_max - 1)^3 must fit into an int. */
  135. #define num_sample_max (1 << ((sizeof(int) * 8 - 1) / 3))
  136.  
  137. /* Table of f(i) = 256 * sqrt(1 + (i/64)^2). */
  138. /* This is good to within 1% or better. */
  139. #define sqrt_index_shift 6        /* scaling of index */
  140. #define sqrt_value_shift 8        /* scaling of value */
  141. private int scaled_sqrt_tab[65] =
  142.    {    256, 256, 256, 256, 256, 256, 257, 257,
  143.     257, 258, 259, 259, 260, 261, 262, 262,
  144.     263, 264, 265, 267, 268, 269, 270, 272,
  145.     273, 274, 276, 277, 279, 281, 282, 284,
  146.     286, 288, 289, 291, 293, 295, 297, 299,
  147.     301, 304, 306, 308, 310, 312, 315, 317,
  148.     320, 322, 324, 327, 329, 332, 334, 337,
  149.     340, 342, 345, 348, 350, 353, 356, 359,
  150.     362
  151.    };
  152.  
  153. private int
  154. flatten_recur(gx_path *ppath,
  155.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3,
  156.   fixed scaled_flat)
  157. {    fixed
  158.       x0 = ppath->position.x,
  159.       y0 = ppath->position.y;
  160. top:
  161. #ifdef DEBUG
  162. if ( gs_debug['2'] )
  163.     dprintf4("[2]x0=%f y0=%f x1=%f y1=%f\n",
  164.          fixed2float(x0), fixed2float(y0),
  165.          fixed2float(x1), fixed2float(y1)),
  166.     dprintf4("   x2=%f y2=%f x3=%f y3=%f\n",
  167.          fixed2float(x2), fixed2float(y2),
  168.          fixed2float(x3), fixed2float(y3));
  169. #endif
  170.     /*
  171.      * Compute the maximum distance of the curve from
  172.      * the line (x0,y0)->(x3,y3).  We do this conservatively
  173.      * by observing that the curve is enclosed by the
  174.      * quadrilateral of its control points, so we simply
  175.      * compute the distances of (x1,y1) and (x2,y2)
  176.      * from the line.  Letting dx = x3-x0 and dy = y3-y0,
  177.      * the distance of (xp,yp) from the line is
  178.      * abs(N)/sqrt(D), where N = dy*(xp-x0)-dx*(yp-y0) and
  179.      * D = dx*dx+dy*dy; hence we want to test abs(N) <= sqrt(D)*F,
  180.      * where F is the flatness parameter from the graphics state.
  181.      * We can do this more efficiently by letting t=dy/dx, and
  182.      * testing abs(N1) <= sqrt(D1)*f, where N1=t*(xp-x0)-(yp-y0) and
  183.      * D1 = 1+t*t.  If dx < dy, we swap x and y for this
  184.      * computation.  This guarantees abs(t) <= 1, which allows us to
  185.      * compute sqrt(1+t*t) by table lookup on the high bits of abs(t).
  186.      *
  187.      * To avoid replacing shallow curves by long straight lines,
  188.      * we also require that abs(dx) and abs(dy) be reasonably small.
  189.      *
  190.      * Note that if scaled_flat is very small, we don't do any of this;
  191.      * instead, we just check whether abs(dx) and abs(dy) are
  192.      * small enough to switch over to sampling rather than dividing.
  193.       */
  194.      { fixed dx3 = x3 - x0;
  195.        fixed adx3 = any_abs(dx3);
  196.        fixed dy3 = y3 - y0;
  197.        fixed ady3 = any_abs(dy3);
  198.        /* We have to be quite careful to ensure that */
  199.        /* none of the multiplications will overflow. */
  200. #define short_max 0x7ff0L
  201. #define reduce_3(ad3, maxv)\
  202.   while ( ad3 > maxv )\
  203.     adx3 >>= 1, ady3 >>= 1,\
  204.     dx3 = arith_rshift_1(dx3), dy3 = arith_rshift_1(dy3)
  205. #define reduce_d(d)\
  206.   for ( shift = 0; (d < 0 ? d < -short_max : d > short_max); shift++ )\
  207.     d = arith_rshift_1(d)
  208.        if ( adx3 > ady3 )
  209.         {    fixed d, dx, dy, dist;
  210.         int shift;
  211.         if ( scaled_flat == 0 )
  212.         { if ( adx3 < int2fixed(num_sample_max / 2) )
  213.           { int n = fixed2int_var_rounded(adx3) * 2 + 1;
  214.             return flatten_sample(ppath, max(n, 3), x1, y1,
  215.                       x2, y2, x3, y3);
  216.           }
  217.           else goto sub;
  218.         }
  219.         else if ( adx3 > scaled_flat << 1 )
  220.           goto sub;
  221.         reduce_3(ady3, short_max);
  222.         d = (scaled_sqrt_tab[(ady3 << sqrt_index_shift) / adx3] * scaled_flat) >> sqrt_value_shift;
  223.         dx = x1 - x0, dy = y1 - y0;
  224.         reduce_d(dx);
  225.         if ( ((dist = ((dx * dy3 / dx3) << shift) - dy) < 0 ?
  226.               -dist : dist) > d )
  227.           goto sub;    /* not flat enough */
  228.         dx = x2 - x0, dy = y2 - y0;
  229.         reduce_d(dx);
  230.         if ( ((dist = ((dx * dy3 / dx3) << shift) - dy) < 0 ?
  231.               -dist : dist) > d )
  232.           goto sub;    /* not flat enough */
  233.         }
  234.        else if ( ady3 != 0 )
  235.         {    fixed d, dy, dx, dist;
  236.         int shift;
  237.         if ( scaled_flat == 0 )
  238.         { if ( ady3 < int2fixed(num_sample_max / 2) )
  239.           { int n = fixed2int_var_rounded(ady3) * 2 + 1;
  240.             return flatten_sample(ppath, max(n, 3), x1, y1,
  241.                       x2, y2, x3, y3);
  242.           }
  243.           else goto sub;
  244.         }
  245.         else if ( ady3 > scaled_flat << 1 )
  246.           goto sub;
  247.         reduce_3(adx3, short_max);
  248.         d = (scaled_sqrt_tab[(adx3 << sqrt_index_shift) / ady3] * scaled_flat) >> sqrt_value_shift;
  249.         dy = y1 - y0, dx = x1 - x0;
  250.         reduce_d(dy);
  251.         if ( ((dist = ((dy * dx3 / dy3) << shift) - dx) < 0 ?
  252.               -dist : dist) > d )
  253.           goto sub;    /* not flat enough */
  254.         dy = y2 - y0, dx = x2 - x0;
  255.         reduce_d(dy);
  256.         if ( ((dist = ((dy * dx3 / dy3) << shift) - dx) < 0 ?
  257.               -dist : dist) > d )
  258.           goto sub;    /* not flat enough */
  259.         }
  260.        else                /* adx3 = ady3 = 0 */
  261.         {    /* (x0,y0) is the same point as (x3,y3). */
  262.         /* This is an anomalous case.  If the entire curve */
  263.         /* is a single point, stop now, otherwise subdivide. */
  264.         if ( x1 != x0 || y1 != y0 || x2 != x0 || y2 != y0 )
  265.           goto sub;
  266.         }
  267.      }
  268.     /* Curve is flat enough.  Add a line and exit. */
  269. #ifdef DEBUG
  270. if ( gs_debug['2'] )
  271.     dprintf2("[2]\t*** x=%f, y=%f ***\n",
  272.          fixed2float(x3), fixed2float(y3));
  273. #endif
  274.     return gx_path_add_line(ppath, x3, y3);
  275.  
  276.     /* Curve isn't flat enough.  Break into two pieces and recur. */
  277.     /* Algorithm is from "The Beta2-split: A special case of the */
  278.     /* Beta-spline Curve and Surface Representation," B. A. Barsky */
  279.     /* and A. D. DeRose, IEEE, 1985, courtesy of Crispin Goswell. */
  280. sub:
  281.     /* We have to define midpoint carefully to avoid overflow. */
  282.     /* (If it overflows, something really pathological is going on, */
  283.     /* but we could get infinite recursion that way....) */
  284. #define midpoint(a,b)\
  285.   (arith_rshift_1(a) + arith_rshift_1(b) + ((a) & (b) & 1))
  286.    {    fixed x01 = midpoint(x0, x1), y01 = midpoint(y0, y1);
  287.     fixed x12 = midpoint(x1, x2), y12 = midpoint(y1, y2);
  288.     fixed x02 = midpoint(x01, x12), y02 = midpoint(y01, y12);
  289.     int code;
  290.     /* Update x/y1, x/y2, and x/y0 now for the second half. */
  291.     x2 = midpoint(x2, x3), y2 = midpoint(y2, y3);
  292.     x1 = midpoint(x12, x2), y1 = midpoint(y12, y2);
  293.     code = flatten_recur(ppath, x01, y01, x02, y02,
  294.         (x0 = midpoint(x02, x1)), (y0 = midpoint(y02, y1)),
  295.         scaled_flat);
  296.     if ( code < 0 ) return code;
  297.    }    goto top;
  298. }
  299.  
  300. /* Flatten a segment of the path by repeated sampling. */
  301. /* n is the number of points to sample, including the endpoints; */
  302. /* we require n >= 3. */
  303. private int
  304. flatten_sample(gx_path *ppath, int n,
  305.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3)
  306. {    const fixed
  307.         x0 = ppath->position.x,
  308.         y0 = ppath->position.y;
  309.     /* We spell out some multiplies by 3, */
  310.     /* for the benefit of compilers that don't optimize this. */
  311.     const fixed
  312.         x01 = x1 - x0,
  313.         cx = (x01 << 1) + x01,        /* 3*(x1-x0) */
  314.         x12 = x2 - x1,
  315.         bx = (x12 << 1) + x12 - cx,    /* 3*(x2-2*x1+x0) */
  316.         ax = x3 - bx - cx - x0;        /* x3-3*x2+3*x1-x0 */
  317.     const fixed
  318.         y01 = y1 - y0,
  319.         cy = (y01 << 1) + y01,
  320.         y12 = y2 - y1,
  321.         by = (y12 << 1) + y12 - cy,
  322.         ay = y3 - by - cy - y0;
  323.     const int
  324.         n1 = n - 1,
  325.         n12 = n1 * n1,
  326.         n13 = n12 * n1;
  327.     fixed ptx = x0, pty = y0;
  328.     fixed x, y;
  329.     /*
  330.      * If all the coefficients lie between min_fast and max_fast,
  331.      * we can do everything in fixed point.  In this case we compute
  332.      * successive values by finite differences, using the formulas:
  333.         x(t) =
  334.           a*t^3 + b*t^2 + c*t + d =>
  335.         dx(t) = x(t+e)-x(t) =
  336.           a*(3*t^2*e + 3*t*e^2 + e^3) + b*(2*t*e + e^2) + c*e =
  337.           (3*a*e)*t^2 + (3*a*e^2 + 2*b*e)*t + (a*e^3 + b*e^2 + c*e) =>
  338.         d2x(t) = dx(t+e)-dx(t) =
  339.           (3*a*e)*(2*t*e + e^2) + (3*a*e^2 + 2*b*e)*e =
  340.           (6*a*e^2)*t + (6*a*e^3 + 2*b*e^2) =>
  341.         d3x(t) = d2x(t+e)-d2x(t) =
  342.           6*a*e^3;
  343.         x(0) = d, dx(0) = (a*e^3 + b*e^2 + c*e),
  344.           d2x(0) = 6*a*e^3 + 2*b*e^2;
  345.      * In these formulas, e = 1/n1; of course, there are separate
  346.      * computations for the x and y values.
  347.      */
  348. #define max_fast (max_fixed / 6)
  349. #define min_fast (-max_fast)
  350.     int fast;
  351.     float dt;        /* only if !fast */
  352.     int i;
  353.     /*
  354.      * We do exact rational arithmetic to avoid accumulating error.
  355.      * Each quantity is represented as I+R/n13, where I is an "integer"
  356.      * and the "remainder" R lies in the range 0 <= R < n13.  Note that
  357.      * R may temporarily exceed n13, and hence possibly overflow.
  358.      */
  359.     fixed idx, idy, id2x, id2y, id3x, id3y;        /* I */
  360.     int rx, ry, rdx, rdy, rd2x, rd2y, rd3x, rd3y;
  361.  
  362.     if_debug6('2', "[2]ax=%f bx=%f cx=%f\n   ay=%f by=%f cy=%f\n",
  363.           fixed2float(ax), fixed2float(bx), fixed2float(cx),
  364.           fixed2float(ay), fixed2float(by), fixed2float(cy));
  365. #define in_range(v) (v < max_fast && v > min_fast)
  366.     if ( n1 < num_sample_max &&    /* so n13 fits */
  367.          in_range(ax) && in_range(bx) && in_range(cx) &&
  368.          in_range(ay) && in_range(by) && in_range(cy)
  369.        )
  370.     {    fast = 1;
  371.         x = x0, y = y0;
  372.         rx = ry = 0;
  373.         /* Fast check for n == 3, a common special case */
  374.         /* for small characters. */
  375.         if ( n == 3 )
  376. #define poly2(a,b,c)\
  377.   arith_rshift_1(arith_rshift_1(arith_rshift_1(a) + b) + c)
  378.             idx = poly2(ax, bx, cx),
  379.             idy = poly2(ay, by, cy),
  380.             rdx = rdy = 0;
  381. #undef poly2
  382.         else
  383.         {    fixed bx2 = bx << 1, by2 = by << 1;
  384.             fixed ax6 = ((ax << 1) + ax) << 1,
  385.                   ay6 = ((ay << 1) + ay) << 1;
  386.             fixed qx, qy;
  387. #define adjust_rem(r, q)\
  388.   if ( r < 0 ) r += n13, q--
  389. #define adjust_rem_loop(r, q)\
  390.   while ( r < 0 ) r += n13, q--;\
  391.   while ( r >= n13 ) r -= n13, q++
  392.             /* We can compute all the remainders as ints, */
  393.             /* because we know they are less than n13. */
  394.             /* bx/y terms */
  395.             id2x = bx2 / n12, id2y = by2 / n12;
  396.             rd2x = ((int)bx2 - (int)id2x * n12) * n1,
  397.               rd2y = ((int)by2 - (int)id2y * n12) * n1;
  398.             idx = bx / n12, idy = by / n12;
  399.             rdx = ((int)bx - (int)idx * n12) * n1,
  400.               rdy = ((int)by - (int)idy * n12) * n1;
  401.             /* cx/y terms */
  402.             idx += qx = cx / n1, idy += qy = cy / n1;
  403.             rdx += ((int)cx - (int)qx * n1) * n12,
  404.               rdy += ((int)cy - (int)qy * n1) * n12;
  405.             /* ax/y terms */
  406.             idx += qx = ax / n13, idy += qy = ay / n13;
  407.             rdx += (int)ax - (int)qx * n13,
  408.               rdy += (int)ay - (int)qy * n13;
  409.             id2x += id3x = ax6 / n13, id2y += id3y = ay6 / n13;
  410.             rd2x += rd3x = (int)ax6 - (int)id3x * n13,
  411.               rd2y += rd3y = (int)ay6 - (int)id3y * n13;
  412.             adjust_rem_loop(rdx, idx);
  413.             adjust_rem_loop(rdy, idy);
  414.             adjust_rem_loop(rd2x, id2x);
  415.             adjust_rem_loop(rd2y, id2y);
  416.             adjust_rem(rd3x, id3x);
  417.             adjust_rem(rd3y, id3y);
  418. #undef adjust_rem
  419. #undef adjust_rem_loop
  420.         }
  421.     }
  422.     else
  423.         fast = 0, dt = 1.0 / (float)n1;
  424.     if_debug4('2', "[2]sampling %s n=%d\n[2]x=%g, y=%g\n",
  425.           (fast ? "fast" : "slow"), n,
  426.           fixed2float(x), fixed2float(y));
  427.     for ( i = 1; i < n1; i++ )
  428.     {    int code;
  429.         if ( fast )
  430.         {
  431. #ifdef DEBUG
  432. if ( gs_debug['2'] )
  433.             dprintf4("[2]dx=%f+%d, dy=%f+%d\n",
  434.                  fixed2float(idx), rdx,
  435.                  fixed2float(idy), rdy),
  436.             dprintf4("   d2x=%f+%d, d2y=%f+%d\n",
  437.                  fixed2float(id2x), rd2x,
  438.                  fixed2float(id2y), rd2y),
  439.             dprintf4("   d3x=%f+%d, d3y=%f+%d\n",
  440.                  fixed2float(id3x), rd3x,
  441.                  fixed2float(id3y), rd3y);
  442. #endif
  443. #define accum(i, r, di, dr)\
  444.   if ( (unsigned)(r += dr) >= (unsigned)n13 ) r -= n13, i++;\
  445.   i += di
  446.             accum(x, rx, idx, rdx);
  447.             accum(idx, rdx, id2x, rd2x);
  448.             accum(id2x, rd2x, id3x, rd3x);
  449.             accum(y, ry, idy, rdy);
  450.             accum(idy, rdy, id2y, rd2y);
  451.             accum(id2y, rd2y, id3y, rd3y);
  452. #undef accum
  453.         }
  454.         else
  455.         {    const float t = dt * (float)(i);
  456.             x = ((ax*t + bx)*t + cx)*t + x0;
  457.             y = ((ay*t + by)*t + cy)*t + y0;
  458.         }
  459.         if_debug3('2', "[2]%s x=%g, y=%g\n",
  460.               (((x ^ ptx) | (y ^ pty)) & float2fixed(-0.5) ?
  461.                "add" : "skip"),
  462.               fixed2float(x), fixed2float(y));
  463.         /* Skip very short segments */
  464.         if ( ((x ^ ptx) | (y ^ pty)) & float2fixed(-0.5) )
  465.         {    if ( (code = gx_path_add_line(ppath, x, y)) < 0 )
  466.                 return code;
  467.             ptx = x, pty = y;
  468.         }
  469.     }
  470.     if_debug2('2', "[2]last x=%g, y=%g\n",
  471.           fixed2float(x3), fixed2float(y3));
  472.     return gx_path_add_line(ppath, x3, y3);
  473. }
  474.  
  475. /*
  476.  *    The rest of this file is an analysis that will eventually
  477.  *    allow us to rasterize curves on the fly, by finding points
  478.  *    where Y reaches a local maximum or minimum, which allows us to
  479.  *    divide the curve into locally Y-monotonic sections.
  480.  */
  481.  
  482. /*
  483.     Let y(t) = a*t^3 + b*t^2 + c*t + d, 0 <= t <= 1.
  484.     Then dy(t) = 3*a*t^2 + 2*b*t + c.
  485.     y(t) has a local minimum or maximum (or inflection point)
  486.     precisely where dy(t) = 0.  Now the roots of dy(t) are
  487.         ( -2*b +/- sqrt(4*b^2 - 12*a*c) ) / 6*a
  488.        =    ( -b +/- sqrt(b*2 - 3*a*c) ) / 3*a
  489.     (Note that real roots exist iff b^2 >= 3*a*c.)
  490.     We want to know if these lie in the range (0..1).
  491.     (We don't care about the endpoints.)  Call such a root
  492.     a "valid zero."  We proceed as follows:
  493.         If sign(3*a + 2*b + c) ~= sign(c), a valid zero exists.
  494.         If sign(a) = sign(b), no valid zero exists.
  495.     Otherwise, we look for a local extremum of dy(t) by observing
  496.         d2y(t) = 6*a*t + 2*b
  497.     which has a zero only at
  498.         t1 = -b / 3*a
  499.     Now if t1 <= 0 or t1 >= 1, no valid zero exists.  Otherwise,
  500.     we compute
  501.         dy(t1) = c - (b^2 / 3*a)
  502.     Then a valid zero exists (at t1) iff sign(dy(t1)) ~= sign(c).
  503.  */
  504.